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ABSTRACT 


The differentiation matrix for a Daubechies-based wavelet basis will be constructed 
and ‘superconvergence’ will be proven. That is, it will be proven that under the 
assumption of periodic boundary conditions that the differentiation matrix is accurate 
of order 2A/, even though the approximation subspace can represent exactly only 
polynomials up to degree M — 1, where M is the number of vanishing moments of the 
associated wavelet. It will be illustrated that Daubechies-based wavelet methods are 
equivalent to finite difference methods with grid refinement in regions of the domain 
where small-scale structure is present. 
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1 Introduction 


The term differentiation matrix was coined by E. Tadmor in his review on spectral 
methods [1]. The term denotes the transformation between grid point values of a 
function and its approximate first derivative. This matrix is a product of three ma- 
trices. 

The first matrix C is constructed as follows: assume that the point values of a 
function f(x) (where a < x < b) are given at N points Xj for 0 < j < N - 1. 
Thus a vector of numbers f(xj ) is given. From this vector one can reconstruct an 
approximation to the function f(x) for every point x in the interval. This approxima- 
tion (denoted by Pjvf) itself belongs to a finite dimensional space - in pseudospectral 
methods it is the global interpolation polynomial that collocates f{xj) and in finite 
differences or finite elements it is a piecewise polynomial. This transformation be- 
tween f(xj) and P N f, defines the matrix C . Of course this matrix depends on the 
special basis chosen to represent P^f. A good example is the Fourier interpolation 
procedure in which the basis is the set of complex exponentials. 

The second matrix D results from differentiating Pyy/ , and projecting it back to 
the finite dimensional space. Thus D is defined by the linear transformation between 
P N f and Pn^PnI- 

The last matrix is the inverse of the first matrix. Basically, since we are given the 
approximation P yv ^Pn f we can read it at the grid points Xj> Thus the differentiation 
matrix T> can be represented as V = C~ l DC. 

In this paper the wavelet differentiation matrix will be examined. As with other 
basis sets, as outlined above, it is a product of three matrices. Under the assumption 
of periodicity of f(x ), however, the matrices C and D commute allowing D to operate 
directly on the vector of numbers f{x 3 ). That is, the differentiation matrix T> is simply 
D : 'D = D, Furthermore, the matrix D differentiates samples of polynomials exactly, 
i.e., the action of D is equivalent to a finite difference operator with order of accuracy 
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depending on the order of the wavelet chosen. 


More precisely, the following outlines the proof of this assertion: 

• Given a periodic function /( x), let C be the mapping from evenly-spaced sam- 
ples of /(x) to the approximate scaling function coefficients on the finest scale: 
C : / — ► (To- Due to the periodicity of /(x), C is circulant in form. 

• Let D be the mapping from the exact scaling function coefficients of f(x) to 
scaling function coefficients of f'(x): D : sq — * so- Once again, due to the 
periodicity of /(x), D is circulant in form. 

• The matrix operator D can differentiate exactly evenly-spaced samples of poly- 
nomials, i.e., D has the effect of a finite-difference operator. The order of exact 
differentiation depends on the order of the wavelet used. 

• All circulant matrices with the same dimensions commute, therefore the oper- 

—f 

ator D can be applied directly to / : 

/' = C-'DCf, 


or simply, 

/' = Df, 

and this will complete the proof. 


This paper contains five sections: 

§1) This introduction. 

§2) Standard definitions of wavelets and scaling functions are given. 

§3) The general approximation properties of wavelets will be discussed along with 
the quadrature formula needed to approximate the scaling function coefficients of a 
function. 
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§4) This is the most important section of this paper. It will be proved that the 
action of D is equivalent to a finite difference operator. 

§5) The results of sections (3) and (4) are combined for the desired conclusion. 

In addition, the following two related topics are explored in the first two appen- 
dices: 

Appendix A) For wavelets supported on (0,3M) it will be shown that / <j>(x)x m dx = 
(f <f>(x)xdx) m 

Appendix B) The moments of the scaling function <f>(x) will be calculated. 
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2 Wavelet Definitions and Relations 


The term wavelet is used to describe a spatially localized function. ‘Localized means 
that the wavelet has compact support or that the wavelet almost has compact sup- 
port in the sense that outside of some interval the amplitude of the wavelet decays 
exponentially. We will consider only wavelets that have compact support and that 
are of the type defined by Daubechies [2] which are supported on [0,2Af 1], where 

M is the number of vanishing moments defined later in this section. 

To define Daubechies wavelets, consider the two functions <f>{x) and ip(x) which 
are solutions to the following equations: 


L - 1 


<f>(x) = \/2 ^2 hk<f>(2x — k ), 

k = 0 

(1) 

L— 1 

V>( x ) = \/2 gk<f>{2z - &)> 

k-0 

(2) 

where <f>(x) is normalized, 

too 

/ <f>(x)dx = 1. 

J — oo 

(3) 

Let, 


<f>{(x ) — 2~* <f>(2~* x — k ), 

(4) 

and 


2T 

I 

H 

1 

1 

CS 

II 

(5) 


where j, k £ Z, denote the dilations and translations of the scaling function and the 
wavelet. 

The coefficients H = {M*=o and G = jt=o are related by g k = {-\) k h L _ k for 
k = 0, L — 1. Furthermore, H' and G are chosen so that dilations and translations 
of the wavelet, x ), form an orthonormal basis of L 2 (R) and so that tp(x) has M 

vanishing moments. In other words, ij>l(x) will satisfy 

SkiSjm = [ il>l(x)rf>F(x)dx, (6) 

J — oo 
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where S k i is the Kronecker delta function. Also, ip(x) = satisfies 

/ oo 

Ip(x)x m dx = 0 , ( 7 ) 

-oo 

for m = 0, ...,M — 1. Under the conditions of the previous two equations, for any 
function f(x) 6 L 2 (R) there exists a set {dj k } such that 


f( x ) = E E d iki/>i( x ), (8) 

j€ZKeZ 

where 

/ oo 

f(x)^{(x)dx. ( 9 ) 

-OO 

The number of vanishing moments of the wavelet ^>(x) defines the accuracy of 
approximation. The two sets of coefficients H and G are known in signal processing 
literature as quadrature mirror filters [3]. For Daubechies wavelets the number of 
coefficients in H and G, or the length of the filters H and G, denoted by L, is related 
to the number of vanishing moments M by 2 M = L. For example, the famous Haar 
wavelet is found by defining H as h Q = h x = 1. For this filter, H, the solution to 

the dilation equation (1), < is the box function: <j>(x) = 1 for x € [0,1] and 

<i>{x) = 0 otherwise. The Haar function is very useful as a learning tool, but it 
is not very useful as a basis function on which to expand another function for the 
important reason that it is not differentiable. The coefficients, //, needed to define 
compactly supported wavelets with a higher degree of regularity can be found in [2]. 
As is expected, the regularity increases with the support of the wavelet. The usual 
notation to denote a Daubechies wavelet defined by coefficients H of length L is Di. 

It is usual to let the spaces spanned by <f>j.(x) and ip k (x) over the parameter k, 
with j fixed, to be denoted by Vj and Wj respectively: 


^• = 


span 

k£Z 


<f>k ( X )> 


span 

kez 


i i ( x )• 


( 10 ) 

(ii) 
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The spaces Vj and Wj are related by [2] 


... CViCVoC V-x C ..., 


( 12 ) 


and that 


Vj = Vj+i©Wj+,- (13) 

The previously stated condition that the wavelets form an orthonormal basis of L 2 (R) 
can now be written as, 

L 2 (R) = ®Wj. (14) 

jez 

Two final properties of the spaces Vj are that 


n^ = {o>. < i5 > 

jez 


and 

u Vj = L 2 (R). ( 16 ) 

jez 

Properties of the Semi-Discrete Fourier Transform (SDFT) of the filter H will also 
be needed. The following definition is not exactly the SDFT but a constant times the 
SDFT: 

(17) 

At=0 

This DFT satisfies the following equation, see [4]: 

\m\ 2 +m +*)\ 2 = i. as) 


Solutions of equation (18) have the following properties, see [2]: 

H{t) = (|(1 + e’ 4 )) M C?(e ,f ), (19) 

where M is the number of vanishing moments of the wavelet and Q is a trigonometric 
polynomial such that, 


|Q(e*| 2 = P( sin 2 (£/2)) + sin 2M (^/2)fi(^ cos £), (20) 
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where 


P(y) = ‘ e " ( M ~i +k ) /, 

and is an odd polynomial such that, 

0 < P{y) + y M R(l/2 - y) 

for 0 < y < 1, and 

SU P (P(y) + y M R( 1/2 - y )) < 2 2(m - 1 >, (23) 

0<y<l v 7 

if A/ > 2 or 

~ - R{x) - t+W\' (24) 

for |x| < 1/2, if M = 1. The important point here is that has a zero of order 

M at £ = 7r. 

Of course, infinite sums and unions are meaningless when one begins to implement 
a wavelet expansion on a computer. In some way one must limit range of the scale 
parameter j and the location parameter k. Consider first the scale parameter j . As 
stated above, the wavelet expansion is complete: L 2 (R) = ® j€Z Wj. Therefore, any 
f{x) € L 2 (R) can be written as, 

jez kez 

where due to orthonormality of the wavelets d{ = f(x)rp j k (x). In this expan- 
sion, functions with arbitrarily small-scale structures can be represented. In practice, 
however, there is a limit to how small the smallest structure can be. This would 
depend, for example, on how fine the grid is in a numerical computation scenario or 
perhaps what the sampling frequency is in a signal processing scenario. Therefore, 
on a computer an expansion would take place in a space such as 

Vo = Wi 0 W 2 ® . . . © Wj © Vj, (25) 


( 21 ) 

( 22 ) 
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and would appear as, 


Pv 0 f(x) = £ s J k <f> J k {x) + £ £ <W( X ), (26) 

k£Z 3=1 keZ 

where again due to orthonormality of the basis functions d 3 k = /_ ^ f{x )V>fc(x), and 
s k = f—oo /( x )^fc( x )- this expansion, scale j = 0 is arbitrarily chosen as the finest 
scale that is needed, and scale J would be the scale at which a kind of local average, 
4> J k (x), provides sufficient large scale information. In language that is likely to appeal 
to the electrical engineer it can be said that <f> J k (x) represents the direct current portion 
of a signal and that represents the alternating current portion of a signal at, 

very roughly, frequency j. As stated above, one must also limit the range of the 
location parameter k. In this paper this is done by assuming that f(x) is a periodic 
function. The periodicity of f(x) induces periodicity on all wavelet coefficients, s k 
and 

This completes the definition of wavelets. The next section will discuss function 
approximation in a wavelet basis. 
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3 Approximating in a Wavelet Basis 


Scaling functions and wavelets were defined in the previous section. The goal of this 
section is to find the coefficients in a wavelet expansion. More precisely, the scal- 
ing function coefficients at the finest scale, So, will be approximated. The key to 
this approximation is the matrix C which maps evenly-spaced samples of a periodic 
function to the approximate scaling function coefficients. This matrix C has the de- 
sirable property of being circulant in form with the ramification that C will commute 
with any other circulant matrix, particularly the derivative matrix RP, the subject of 
section (4). An example of C is given at the end of section (3.2). 

The scenario for this section is as follows: let the finest scale be scale j = 0, i.e., 
at this scale all relevant small scale structures in the function have been captured and 
represented. One seeks an expansion of a function f(x) in terms of <f>° in the space 
Vo- With the projection Py 0 defined as Py 0 : L 2 (R ) — ► Vo such an expansion has the 
following form: 

A './(*) = £ (*), (27) 

kez 

where due to the orthonormality of the basis functions, <f> ( -(x)(f> < j(x)dx, the 

coefficients s° are given by 

/ OO 

f(x)<f>° k (x)dx. (28) 

-OO 

Once the s k have been found one would usually then find the scaling function and 
wavelet function coefficients at more coarse scales. This can be done by using equation 
(29) to get s{ for j = 1 ,..., J and by using equation (30) to get d J k for j = 1 ,..., J. 
These equations are derived respectively from equations (1) and (2), see [4], [5], 

2 M 

S k ~ XI ^" S n+2fc-25 (29) 

n=l 

and 

2 M 

dk ~ XI 9nS J n+ 2k-2- (30) 
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In this section, however, the decomposition onto more coarse scales will not be cal- 
culated. The important step for this section is the approximation of the integral 
s° k = f{x)(fP k {x)dx. Let er° k denote the approximation to s° k . The quadrature for- 
mula for this integral encompasses the approximation properties of scaling functions, 
and hence wavelets. 

This section contains 3 subsections: 

§3.1 ) The approximation properties of scaling functions will be discussed and the 
quadrature formula to estimate the integral f(x)<f>° k (x)dx will be derived. 

§3.2 ) An example using the results from section (3.1) is given for the Daubechies 
wavelet D&. 

§3.3 ) The example from section (3.2) leads to a circulant matrix for the matrix 
C. Circulant matrices will be defined and the ramifications of circularity will be 
discussed. 

3.1 Quadrature Formula for Scaling Function 

In this subsection the coefficients s° k will be approximated. Before stating the appro- 
priate quadrature formula, however, the order of accuracy of a wavelet approximation 
is discussed. 

3.1.1 Approximation Properties of Scaling Functions 

This subsection comes essentially from [6]. The approximation properties of scaling 
functions are determined by the Discrete Fourier Transform of the filter H. That is, 
if 

= (31) 
** k= 0 

has a zero of order M at £ = ir then there are a number of consequences: 

1. The polynomials 1, x, ..., are linear combinations of the translates of the 

scaling function (j>° k . 
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2. Smooth functions can be approximated with error 0(h M ), where h denotes the 
grid size, i.e., there exist a set s k , where j is fixed, and there exists a constant 
C such that 

I /(*) - E*Wf(*)l < CA"|/<">( I )|, (32) 

k 

where the norm | • | is the norm. 

3. The associated wavelet has M vanishing moments, 

J x m ip(x)dx = 0 

for m = 0 , M — 1. 

Other ramifications can be found in [6]. These approximation properties determine 
the accuracy of the quadrature formula used to approximate the scaling function 
coefficients sq which is derived in the following section. 


3.1.2 Derivation of Quadrature Formula 

It is important to note that all wavelets in this paper are of the usual Daubechies 
type, i.e., the support of a usual Daubechies wavelet D 2 M is [0, 2A/ — 1] where M is 
the number of vanishing moments of the wavelet. For this subsection this support size 
is particularly important to keep in mind because there does exist an orthonormal 
family of wavelets which are supported on [0, 3 M - 1] and which have a very simple 
quadrature formula based on the vanishing moments of the wavelet (see appendix A) 
but this is not the wavelet being used in this paper. 

Given the approximation properties of the scaling function from the previous 
subsection, one can now seek a quadrature formula which is exact when f(x) is a 
polynomial up to order M - 1: f(x) = p(x) G P W -i. That is, there exist a set of 
coefficients {c/}^ 1 such that 


roo M - • ! 

/ p{x)4>° k dx = Y, c iP( l + *), 
1=0 


(33) 
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for p(x) G Pm-i- If the integral is shifted the above equation becomes, 

roo 

I p{y + k)<f>o(y)dy = E cip{l + k). 

J-OO J =0 


( 34 ) 


More simply, the coefficients can be found [9] by solving the following linear 

system: 

[ x m <f>(x)dx = E l m c h (35) 

J -°° 1=0 

for m = 0, 1 , M — 1, and the coefficients {c,}^ 1 provide the desired quadrature 

formula. That is, the coefficients, (To, which approximate so are found from, 

Af-l 

*1 = E «/(* + *)■ ( 36 ) 


When placed in matrix form the coefficients {c/}^ 0 1 yield the circulant matrix C. 
A more thorough discussion of circulant matrices will be given in §(3.3) after the 


example of the next subsection has been completed. 

Note that since the above derived quadrature formula is exact for p(x) € Pm-i 
the coefficients erg approximate the coefficients s° k with error of order M. Also, note 
that the derivation of the quadrature coefficients depends only on the moments of the 
scaling function, x m <j>{x)dx . In the next subsection, the moments of the scaling 

function will first be calculated and then the coefficients will be found for 

the D$ wavelet. The wavelet Z) 6 is chosen for no other reason than that D 2 and D 4 
receive considerable attention from other sources and that D$ is slightly less trivial 
than D 2 and D 4 while remaining manageable. 


3.2 Example with the Daubechies Wavelet Dq 

Recall from the previous subsection that the immediate goal is to approximate the 
scaling function coefficients of & function at scale j = 0. Specifically, in this section 
the objective is to derive the matrix form of the mapping from evenly-spaced samples 
of a periodic function f(x) to the scaling function coefficients on the finest scale s° k . 
The example will be for the Daubechies wavelet Z) 6 . Comparable results for the 
wavelets D 4 and Dg are presented in appendix B. 
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Recall from the previous subsection that in order to calculate the coefficients 
{ c /}/=o t'he moments of the scaling function <f>(x) must first be known. Let Mi be 
the I th moment of the scaling function <f>(x) y 


Mi = J <f>(x)x‘dx, 

and let \i\ be the I th moment of the filter h k , 

(37) 

W = k‘h k . 

k 

The zero-th moment, M 0 , of <f>(x ) is 1 by the normalization of <f>(x ): 

(38) 

!= 

II 

fr 

II 

> — ‘ 

(39) 

The zero-th moment of the coefficients h k is found by integrating the dilation equation 
which defines (f>{x ): 

J <j>(x)dx = Y, h k j - k)dx, 

and let y = 2x — k to get, 

(40) 

1 = \ H h * J <i>(y) d y, 

which implies, 

(41) 

*> = !> = 2. 

(42) 


k 


That is, the zero-th moments Mo and fi o &re the same for all Daubechies wavelets. 
Higher moments for fii can be found by straight-forward calculation using the coeffi- 
cients provided by Daubechies [2]. The higher moments, Af/ for / > 0, for the scaling 
function can be found from the following equation which is derived in appendix B: 

= (ir +1 £ ( 7 ) (43) 

For the current example only the moments M 0 , M u and M 2 are needed: M 0 = 1, 
Mi = an< ^ M 2 = |((/xi) 2 + fi 2 ). Given these three moments the coefficients 

{c,}^ 1 can be found from 

M - 1 

l m Ci = 

1=0 



J x m <f>(x)dx 


(44) 
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i Mi Pi c, 

0 1 2 .1080 

1 .8174 1.6348 .9667 

2 .6681 1.3363 -.0746 


Table 1: Scaling function and filter moments for Daubechies 6 wavelets. 


for m = 0, 1, ..., M - 1. Specifically, for the D 6 wavelet the linear system in matrix 
form is, 

(l I 2 ) ( 2 ) = ( Mi ) , (45) 

V 0 1 4 / V c 2 / \m 2 ) 

which has the solution cq = .1080, c x = .9667, and c 2 = -.0746. In tabular form, the 


complete results for D& are, 

Recall that the quadrature formula used to approximate the scaling function has 


the form, 


<7k 


M - 1 


C '/(* + *)• 

;=o 


If the function f(x) is periodic 


then in matrix notation the above operation is ctq 


(46) 
= Cf 


where C for D$ and on a grid of 6 points is, 


.108 .967 -.075 0 0 0 \ 

0 .108 .967 -.075 0 0 

0 0 .108 .967 -.075 0 

0 0 0 .108 .967 -.075 * 

-.075 0 0 0 .108 .967 

< .967 -.075 0 0 0 .108 / 


(47) 


The important point here is that the above matrix is circulant. The ramifications 
of circularity are very important for the thesis of the paper. The definition of circulant 
matrices and the properties that they are imbued with is the subject of the next 


subsection. 


3.3 Circulant Matrices 

Strang [7] defines a circulant matrix as a constant-diagonal matrix which is ‘periodic, 
since the lower diagonals fold around to appear again as the upper diagonals.’ A 
thorough discussion of circulant matrices is given by Davis [8]. Circulant matrices 
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have the wonderful property that they can all be diagonalized by the same matrix, the 
Fourier matrix: An NxN Fourier matrix has as its ij-th element the entry u,(‘-i)0-i) 
where w N = 1. The most important ramification for this paper is that matrices which 
can be diagonalized by the same matrix commute. That is, the matrix C from the 
previous subsection will commute with the matrix which will be derived in section 
( 1 . 4 ). 

In general, circulant matrices arise whenever one is performing the matrix version 
of periodic discrete convolution. In numerical analysis periodic discrete convolution 
arises whenever one differentiates the evenly-spaced samples of a periodic equation 
which has constant coefficients. Let us be a bit more precise and illustrate how the 
operation of periodic discrete convolution yields a circulant matrix by stating the 
following theorem; 

Theorem: A finite-length filter of length M applied to N evenly-spaced samples 
of periodic function, where N > M, will in matrix form yield a circulant matrix. 

Proof: First of all, let the notation remain as above: co,cj, j will represent 

the finite-length filter and /o, /i, ..., /jv-i will represent the evenly-spaced samples of 
one period of the periodic function f(x). Of course, the samples of f(x) are periodic 
with period N. The application of the filter c on the samples of f(x) is the convolution: 

M-\ 

a k = C lfk-h 
1-0 

where /, is the renaming of the elements of /,* so that the previous convolution is 
the same as the following expression. Furthermore, keep in mind that / t and /, are 
periodic with period N. 

M - 1 

a k — ^2 c ifi+k- 
1=0 

Using the modulus function to keep the indices of / t - within one period, i.e., keep 
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0 < i < N — 1, the above equation can be written as, 

0k = c lfmod(l+k,N) • 

1=0 

Now, shift the indices by letting j = l + k to get, 

*+(M- 1) 

<?k = 53 c j-kfmod(j,N)‘ 

j=k 

If the length -M filter c is now ’padded’ at the end with zeros so that it is now a 
length - N filter then the above equation can be rewritten as, 

AT— 1 

= 2J c mod(j—k,N)fj' 

3 = 0 

This is exactly a matrix multiply B = Cf where the ij — th element of the matrix 
C is Cnody.i'N), and this is the definition of a circulant matrix. This completes the 
proof./ / 

Note that the difference between a circulant matrix and a Toeplitz matrix is the 
wrapping around effect of the diagonals introduced by the use of the modulus function 
for the circulant matrix. That is, a circulant matrix is a special case of a Toeplitz 
matrix where the constant diagonals are periodic. 

Before leaving the discussion on circulant matrices let one more interpretation be 
noted: to say that circulant matrices commute is to simply restate the important re- 
sult from signal analysis that convolutions commute. That is, if one has two sequences 
c and r, which in the current scenario are periodic, then the order of convolution does 
not matter. This is easily proved with the Fourier transform: 

c*~r = cr = rc = r*~c, (48) 

where r denotes the Fourier transform of r and V denotes periodic convolution. For 
this paper, c, of course, would be the quadrature operator and r would be the scaling 
function derivative operator which is the subject of section (4). In matrix notation 
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equation (48) is nothing more than, 


CH° = H°-C (49) 

In this section a quadrature formula has been found to approximate the scaling 
function coefficients of a given function, f(x). In matrix form this quadrature for- 
mula leads to a circulant matrix assuming f(x} is periodic. In the next section the 
wavelet derivative operator will be derived, and it will be shown that, once again, 
the assumption that f(x) is periodic leads to an operator which in matrix form is 
circulant. 
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4 Derivative based on Wavelets 


In the previous section the mapping from evenly-spaced samples of a periodic function, 
f(x), to the approximate scaling function coefficients on the finest scale, cr k , was given. 
Recall that cr£ denotes the approximation to the exact scaling function coefficient s 3 k 
at scale j and position k. The mapping is nothing more than a quadrature formula 
which is exact when f(x ) is equal to a polynomial up to order M — 1, where M is 
the number of vanishing moments of the wavelet. The question now is how does one 
represent the derivative of f(x) in the wavelet basis given that the wavelet expansion 
of f(x) is already given. 

The answer is given in the following subsections: 

§4.1 ) A function f(x) will be expanded in a wavelet basis and the expansion will 
be differentiated. 

§4.2 ) The results from Beylkin [9] on derivative projections will be given. 

§4.3 ) First, it will be noted that one can differentiate a wavelet expansion at any 
level of a wavelet decomposition and achieve the same derivative. Second, explicit 
wavelet decomposition will be performed accompanied by the appropriate differenti- 
ation matrix. 

§4.4 ) The similarity between the wavelet-based derivative coefficients and finite 
difference derivative coefficients will be noted, and it will be shown that when one 
differentiates the wavelet expansion of a periodic function that the effect on the 
original function samples is equal to finite difference differentiation. 

4.1 Expansion in a Wavelet Basis 

The goal now is to find the wavelet and scaling function expansion of a periodic 
function f(x). Given f(x) € L 2 (R) one first projects onto the arbitrarily chosen 
finest scale j = 0 of the scaling function <f>° k (x) which generates the space V 0 , i.e., let 
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( 50 ) 


Pv b be the projection from the space L 2 (R) to the space Vo, Pv a : L 2 (R) — » Vo: 

PvJ(x) = E »M(*), 

k = 0 

where due to the orthonormality of <f>° k over k in Vo, 

/ OO 

f(x)<t>° k (x)dx. (51) 

-OO 

Note that in the introduction the projection denoted by P would be the same as 
Pv 0 using notation that is more amenable to wavelets. The derivative of Pv 0 f{x ) is, 

—PvJ(x) = £ s° k <f° k (x). (52) 

ax k=o 

Of course, ^P Vo f(x) is not in V 0 and must be projected onto Vo. First define the 
inner product < f,g > on L 2 (R) by 

/ OO 

f(x)g(x)dx. ( 53 ) 

-OO 

Now the projection of £P Vo f(x) onto V 0 is, 

Pv. j- x p vj{x) =E< 2 i p vj’ *? > #?(*). ( 54 ) 

or, 

Pv o4~Pv 0 f(x) = s °k < > 4>K X )- ( 55 ) 

aX 1—0 k=0 

/ 

The inner product < > is one of the results provided in [10]. 

In the previous paragraph f(x) was expanded in a scaling-function expansion at 
the finest scale j = 0. Now f(x) will be expanded in terms of scaling functions and 
wavelets at scale j = 1 . Recall that Vo = V\ © W\ . Now one must project from L 2 (R) 
onto V\ and from L 2 (R) onto W\. Let both projections be denoted simultaneously 
by /View*. That is, P Vi(BWl : L 2 (R) Let Pv^Wifix) be the projection 

°f f(x ) on V\ © W\. Therefore, the expansion for Pv^Wif{x) is, 

N/ 2-1 N/ 2-1 

Pv l9 wJ(x)= 5 fc^fc( X )+ d k^k( X )^ (56) 

k=0 k=0 
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where due to the orthonormality of the basis functions <f>\{x) and ipl(x) the coefficients 
aj and d* are given by 

s l= [ f( x )<t>l( x ) dx ’ ( 57 ) 

OO 

and 

d\ = / f(x)xl>l(x)dx. (58) 

J — OO 

The derivative of Pv^Wifix) is 

, N/ 2-1 , N/ 2-1 

3p/Vi©wi/(^) = ■ s *$t( x ) + 13 

ax k= 0 k = 0 

Once again, the derivative of Py, @vv', /(x) does not belong to Vj ® W \ , and must, 
therefore, be projected back onto this space. The projection is, 


PViQWi /( X ) = (®®) 

N/ 2-1 N/ 2-1 

E E 4 < > <f>H x ) 

1=0 k = 0 
N/2—1 N/ 2-1 

+ E E 4 < > 0/ (*) 

;=o Jt=o 
N/ 2-1 N/ 2-1 

+ E E d k < ^ki^l ^ $( X ) 

;=o Jt=o 
N/2—\ N/ 2-1 

+ E E 4 < > #( x )- 

i=0 fc=o 

The four inner products < </[, $ >, < V'/ >» < V’jti 4>] >, an d < V’/ > are 

the key to finding the derivative of a wavelet expansion, and are provided in [9]. An 
outline of the derivation of these inner products is given in the next section. 


4.2 Wavelet Coefficients of the Derivative 


An arbitrary wavelet expansion of a function might contain wavelet coefficients and 
scaling coefficients at many scales. In [9] the projection coefficients that map from 
scaling function coefficients and wavelet function coefficients at a given scale to the 
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derivative scaling function coefficients and wavelet function coefficients at the same 
scale are derived. The matrix elements of these projections are computed from, 



too / 

2 J «i-/ = a\ { = 2 -2j / ip(2~ : ’x — i)^>( 2~*x — l)dx, 

J — oo 

(61) 


too , 

2 J bi-i = = 2 2j / V>( 2~-’x — i)<j>(2~*x — l)dx, 

J — OO 

(62) 


too , 

2 = 4 = 2 _2j / <f>{ 2~ 3 x - i)i/>{ 2~ 3 x - l)dx, 

J — OO 

(63) 


too , 

2 _J r,_/ = 4 = 2~ 2j / <f>(2~ :l x - 2~ 3 x - l)dx. 

J — oo 

(64) 

Since the above projections are always at a fixed scale, j, the projection coefficients 

are simply, 

0 
II 

1 

sT* 

e- 

(65) 


roo A 

bi= 4>(x - l)—(f>(x)dx, 

J — OO dx 

(66) 


too d 

ci = <t>(x - /)— rj>(x)dx, 

J — oo dx 

(67) 


d 

r i — <f>(x - l)—<f>(x)dx, 

J —oo dx 

(68) 

for l E Z. 

Furthermore, using the dilation equation which defines <^(x), <f>(x) = 

J2k hk<f>{2x - 

- A:), and the equation which defines ip(x), i/>(x) = 'iZk 9k<t>(%x 

— k), the 

first three of the above four equations become, 



L- 1 L— 1 

a i — ^2 Y2 9k9l r 2i+lc-l 
k = 0 1=0 

(69) 


L- 1 i,-l 

b%—^2^2 9khtr2i+k-i 

lt=0 1=0 

(70) 


L- 1 L-l 

c t = ^2 hk9i r 2i+k-i- 

k = o ;=o 

(71) 


It is apparent from the above equations that the coefficients r/ contain all the infor- 


mation concerning the derivative. The coefficients r, can be found [9] from solving 
the following system of linear algebraic equations: 

2 i/2 

n = 2 (r 2 / + - ^ oi 2 k—i ( 2 * 2 / — 2 fc + 1 + r 2 i+ 2 k-i)), (72) 

z k = l 
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and 


where 




L— 1 — n 

= 2 ^ ^ hihi+m 

t—O 


( 73 ) 


(74) 


for n = 1, ..., L — 1. The proof of the above proposition can be found in [9]. 

This section has given a brief outline of the derivation of the wavelet derivative 
projection coefficients. It is important to note that all the information for the wavelet 
derivative is contained in the coefficients {r/}, and this point will be explored more 
in next section. 

4.3 Derivative at Scale Zero of Scaling Function Only 

Wavelet derivatives can be calculated at any level of a wavelet decomposition. The 
result will, of course, always be the same. That is, recall the relation from §(2), 
Vj = Vj+ 1 © Wj +1 . As stated before, it is the convention of this paper to let Vo represent 
the finest scale. Using the above relation, one could decompose Vo any number of 
times. One decomposition yields Vo = W\ © Vi, and a second decomposition yields 
V 0 = Wi@W 2 <$ V 2 . One could calculate the wavelet derivative in any one of these 
spaces. Once again, the goal of this paper is to understand the essence of a wavelet 
derivative, and since the derivative is the same regardless of the decomposition of the 
space, one should choose the simplest approach and calculate the derivative at scale 
j = 0 using only the scaling function coefficients. 

This subsection contains four parts: 

1. New notation will be introduced. 

2. Wavelet decompositions and differentiation matrices will be given for the space 
Vo as well as comments on data compression in this space. 

3. Wavelet decompositions and differentiation matrices will be given for the space 
W\ © V\ as well as comments on data compression in this space. 
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4. Wavelet decompositions and differentiation matrices will be given for the space 
Wi © W 2 © W 3 © V3 as well as comments on data compression. 


4.3.1 New Notation 


To simplify the presentation, matrix notation will be used whenever possible. Begin 
by defining the matrix version of equations (29) and (30). Recall that these equations 
are 


n=2Af 
n= 1 


and 


n=2M 


~ 5 r "' S n+2Jk-2- 

n= 1 


Denote the decomposition matrix embodied by these two equations by Ptf+x where 
the matrix subscripts denote the size of the matrix, and the superscripts indicate 
that P is decomposing from scaling function coefficients at scale j to scaling function 
and wavelet function coefficients at scale j + 1. As before, let sj contain the scaling 
function coefficients at scale j. (Note: When vector notation is used the scale is 
given as a subscript, otherwise the location k is the subscript and the scale is the 
superscript.) P therefore maps Sj onto s J+1 and dj+\ : 


pj.j+i 

r NxN 


Sj+1 

dj + 1 J 


(75) 


Note that the vectors at scale j + 1 are half as long as the vectors as scale j. To 
illustrate further, suppose the wavelet being used is the four coefficient D 4 wavelet, 
and suppose one wants to project from 8 scaling function coefficients at scale j to 4 
scaling function coefficients at scale j + 1 and 4 wavelet coefficients at scale j © 1. 
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The decomposition matrix in this case is, 

hi hi h$ h+ 0 0 0 0 

0 0 ht h 2 h 3 h 4 0 0 

0 0 0 0 hi h 2 h 3 h 4 

pjj+1 _ I h 3 h 4 0 0 0 0 hi h 2 

8x8 ' ' Ji ft J3 J4 0 0 0 0 

0 0 gi g 2 03 94 0 0 

0 0 0 0 gi g 2 gz g* 

L 93 94 0 0 0 0 gi g 2 

Other decomposition matrices of different sizes will have the same structure as the 


( 76 ) 


above matrix. 

For a bit more matrix notation, let the four matrices A 3 NxN , B 3 NxN , C 3 n x/V, an( l 
R 3 NxN contain the derivative projection coefficients defined in §(4.2) where, again, 
the subscripts denote the size of the matrix and the superscript denotes the scale on 
which the derivative projection is acting. The elements of the four matrices are, 


A <-» a<j = a,_j, 

B bij = &t_j, 
C <-> Cij — Cj — j , 


and 


R r,j = r, w , 


and the mappings performed by the matrices are, 

A 3 : dj ^ dj, 

B 3 : sj - d } , 

C 3 : di - l h 

R 3 : Sj -► sj, 

where Sj and dj, as before, define the scaling and wavelet coefficients at scale j, and 

-f 

Sj and dj denote the coefficients of the expansion of the derivative of a function which 
is initially defined by an expansion in Sj and dj. 

This concludes the new notation. 
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4.3.2 Wavelet Expansion and Derivative in Vo 

As stated previously, one can calculate the derivative of a wavelet expansion at any 
level in the wavelet decomposition. This subsection will explore the first of three of the 
options. To be explicit, suppose that a periodic function f(x ) has been approximated 
on a grid with 16 scaling function coefficients to get s 0 , and for the current argument 
assume that the coefficients have been calculated exactly, i.e., the notation s 0 will be 
used instead of <f 0 . Furthermore due to the periodicity of f(x) the coefficients s 0 will 
also be periodic. The coefficients of the expansion of ~^f{x) in Vo are found from sq 
by an application of the previously defined matrix f?? 6xl6 : 



which completely defines the derivative of f(x) in the scaling function basis at scale 

j = o. 

For data compression purposes, the space V 0 is not a good space to work in. That 
is, the coefficients s*o represent the equivalent of a local averages. In a wavelet basis, 
it is often true that the coefficients of local high-frequency oscillations are small and 
can be set to zero without altering the character of the function being represented, 
but the coefficients of local averages usually represent essential information. 
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4.3.3 Wavelet Expansion and Derivative in W\ ® Vi 

Consider now a decomposition of the vector of scaling function coefficients s 0 onto 
the scaling function and wavelet coefficients at scale j = 1 by an application of the 
matrix Pw xl6 : 



As before, we have 16 basis functions in our space which is now Vi ® W\ rather than 
Vo- In order to calculate the coefficients of the derivative expansion in V\ © W\ the 
following projections are calculated: 

= #8x8 ‘ si + Cl xS • d u (^9) 


= ^8x8 * ^1 + #8X8 ' Si, 


where A , B, C, and R were all defined in the previous subsection. A more concise 
way to represent the derivative projections is in matrix notation: 


Sl _ #8x8 ^8x8 . 

d\ . #8x8 -^8x8 . . ^1 


(81) 


If one now applies the matrix (P X q\ ie) T (T denotes transpose and hence inverse for 
this unitary matrix) to the derivative coefficients at scale j = 1 then one gets the 
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derivative coefficients at scale j = 0: 

io]=(C*.«) T - ? , ( 82 ) 

L ** 1 

and one gets exactly the same coefficients as before when the matrix R j 6xl6 was 

applied to s 0 . To emphasize, the derivative calculated at scale j = 0 and the derivative 

calculated at j = 1 yield exactly the same result. The importance of this observation 

is that in order to understand the essence of the wavelet derivative one need only be 

concerned with the action of the matrix Rf^ X ff on the vector So- 

For data compression the space W\ © V\ is a fair space to work in. The coefficients 

s\ of the basis functions in V\ represent local averages just as the coefficients of 

the basis functions in the space Vo do. However, the basis functions in Vj have 

broader support than the basis functions in Vo and therefore represent averages over 

a larger amount of data (twice as much data to be exact). Therefore, once again 

— # 

the coefficients s*x usually carry essential information. The coefficients d\ of the basis 
functions in the space W\, on the other hand, carry information concerning local 
oscillations. That is, if the function being represented, f(x ), is globally smooth then 
the coefficients d\ will be near zero and can be set exactly to zero without altering the 
character of f(x). In the solution of nonlinear partial differential equations where a 
sharp gradient, or shock, can develop, the coefficients d\ away from the shock would 
be close to zero whereas the coefficients near the shock would be large. Therefore, 
representing a function in W\ © Vi is more versatile than simply staying in the space 
Vo. Versatility continues to be enhanced as one decomposes into more and more 
wavelet subspaces as in the next and final scenario. 

4.3.4 Wavelet Expansion and Derivative in W\ © W 2 © W 3 © V 3 

Up to now our basis functions have all been at the same scale, i.e., initially our 
basis functions were contained in Vo, and in the second scenario the basis functions 
were contained in Vj and W\. In this subsection, however, the basis functions will 
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be contained in spaces at three different scales: W\, W2, W3, and V3. The full set 
of coefficients in this case and all the appropriate decompositions leading to these 
coefficients are, 



In matrix form the projection onto the coefficients of the derivative of the expansion 
is, where the matrix will be labeled M, 
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and M performs the following mapping: 


M 


„3 1 


\ 4 ] 

4 

3 


4 

4 


r u 1 
J3 

'4' 

MS 

[-<*?- 


«1 

4 1 

4 


4 

4 



1.4 J 


[4 J 

[41 


\d\] 

4 


4 

4 


4 

4 


4 

4 


4 

4 


d\ 

4 

[dl\ 


4 

.4 J 


( 85 ) 


For data compression, this is the most useful set of subspaces. The space now is 
represented as W\ 0 W 2 ® W$ ® V3. For the same reasons as before the coefficients 
of basis functions in the subspace V3 cannot be ignored. It is likely, however, that 
the function f(x) being represented is smooth in most of the domain allowing one 
to disregard the majority of the coefficients of the basis functions in the subspace 
Wi ® W 2 ® W3. In fact, it is more likely that the coefficients for the basis functions 
in W\ will be negligible than for the coefficients for the basis functions in W3. This 
is because the basis functions in W3 have larger support than the basis functions in 
W 2 and Wi. 


In summary, an attempt has been made to illustrate that the derivative coefficients 
of a scaling and wavelet expansion can be calculated at any scale. The proper set of 
wavelet subspaces depends on the problem at hand. The goal for this author is to 
understand exactly what wavelets are and what they are doing, therefore, scale j = 0, 
i.e., the space Vo, provides the clearest scenario in which to work without sacrificing 
essential properties of wavelets. 

Given, now, that it is sufficient to work on scale j = 0 to understand exactly 
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what the wavelet derivative does, one must understand the ramifications of applying 
the matrix BP to the vector so. In the next subsection the similarity between the 
above defined matrix BP and finite difference formulas for taking the derivative will 
be explored. 

4.4 Wavelet Derivatives and Finite Difference 

As the previous subsection illustrated, the essential properties of the wavelet deriva- 
tive are contained in the elements of the matrix B. Recall that B is the matrix form 
of the mapping from so to s'o- The surprising property that the matrix B exhibits is, 
however, that it can also differentiate evenly-spaced samples of a function. That is, 
B acts as a finite-difference operator when applied to the samples of a function. 

This subsection is in three parts: 

1. The similarity between wavelet derivative coefficients and finite difference coef- 
ficients is noted. 

2. The finite difference accuracy of the coefficients {r/} derived in [9] will be illus- 
trated, and it will be proved in general that the coefficients {r/} can differentiate 
polynomials exactly up through order 2 M for coefficients {r/} that were derived 
for Daubechies wavelets D 2 m- 

3. In the finite element method under certain conditions one achieves a very high 
order of accuracy called ’superconvergence.’ In wavelet differentiation a similar 
phenomenon is encountered. This phenomenon is defined and a short explana- 
tion is offered. 

4.4.1 Finite Difference Coefficients 

First of all, it is useful to simply note the similarity between the coefficients of centered 
finite difference formulas and the coefficients used to construct the matrix B. The 
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Table 3: Scaling function derivative convolution coefficients for Daubechies wavelets. 

following is a table of centered finite difference coefficients and the order of accuracy 
of the approximation to the derivative: 

Recall that the elements of the matrix R derived in [9] provide the transformation 
from scaling function coefficients of a function to the scaling function coefficients of 
the derivative of the same function. The elements of R for the D<i and D\ wavelet 
derivatives are, as is shown in the following table, exactly the same a s the coefficients 
for the 2-nd and 4-th order centered finite difference formulas. Note that the wavelet 
filters become quite long with increasing order. Therefore, only the right side of the 
filter will be shown keeping in mind that these filters are antisymmetric: 

The fractions for the Dg and Dg wavelets are exact but complicated and provide 
little insight. Compare the following decimal representations of the 6-th and 8-th 
order finite difference operators to the decimal representations of Dg and Dg. Once 
again, only the right-hand side of these antisymmetric filters is shown: 

The coefficients for the Dg and Dg derivatives are not the same as the coefficients 
for the 6-th order and 8-th order centered finite difference derivatives, but the dif- 
ferences are not large. Surprisingly, however, Dg has the same accuracy as the 6-th 
order finite difference operator, and Dg has the same accuracy as the 8-th order finite 
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FD-6 and D 6 

Coefficients 

FD-6 

0 .750 -.150 .017 

D 6 

0 .745 -.145 .015 .0003 


Table 4: Comparison between numerical values of optimal 6th order centered finite 
difference coefficients and Daubechies 6 scaling function derivative convolution coef- 
ficients. 


FD-8 and D$ Coefficients 

FfT8 0 .80 -.20 .038 -.0036 

P 8 0 ,79 -.19 .034 -.0022 -.0002 .0000008 

Table 5: Comparison between numerical values of optimal 8th order centered finite 
difference coefficients and Daubechies 8 scaling function derivative convolution coef- 
ficients. 

difference operator. §4.4.2 will establish this accuracy. 


4.4.2 Finite Difference Accuracy 

To establish the finite-difference accuracy of the wavelet-based differentiation coeffi- 
cients note that a centered-finite-difference derivative approximation with 2 K anti- 
symmetric, r* = — r_jt implying r 0 = 0, coefficients, can be written 

f(xj) ~ £ r* {fj+k - fj-k )• (86) 

fc=i 

If the above equation is exact for f(x) = x q for q = 0, ..., N but not for q — N + 1 
then the equation is said to be N-th order accurate. Therefore, one must check to see 
if 

= £ rk ( x< j+k - ( 87 ) 

k = i 

when f(x) = x q . To simplify, one can let Xj = j and check the following: 

r *(0? + k ) q ~ ti ~ k ) 9 )- ( 88 ) 

k = 1 

Now, treating the coefficients derived in [9] as nothing more than finite- difference 
coefficients one can check the accuracy. The following table contains the results of 
applying the coefficients from [9] to various polynomials: 
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Wavelet Derivative 

Exact up to 

But not for 

d 2 

X 2 

x 3 

D 4 

x* 

X s 

D 6 

X 6 

x 7 

D s 

X 8 

X 9 

Djo 

x 10 

X 11 

D 12 

X 12 

X 13 

D 14 

X 14 

X 15 

D 16 

x 16 

x 17 

D \ 8 

X 18 

X 19 


Table 6 : The degree of polynomials differentiated exactly by various Daubechies 
scaling function derivative coefficients. 

The pattern in this table is obvious and leads to the following theorem: 
Theorem: If <j>{x) is the scaling function for the Daubechies wavelet denoted 
by D 2M , where M is the number of vanishing moments of the wavelet, then the 
coefficients {r;} derived from 

/ oo d 

<f>(x — l)—<j>(x)dx 
-oo dx 

and applied the to evenly-spaced samples of a function act as a finite difference 
derivative operator of order 2 A/. 

The proof of this theorem requires two results, as well as the Fourier Transform 
of <j>(x). First <^(£) will be found followed by the two results which are needed and 
stated in theorem form. Perhaps the first proof would suffice assuming the second 
result is well-known. However, to be complete the second result is also proved. 

Fourier Transform of <f>(x ): Recall, 

L - 1 

<j>{x) = s/2 ^2 hk<f>(2x — k). 
k = o 

Define the Fourier Transform of <f>(x) as, 

A t 00 

m = / <t>(x)e' ix dx. 

J — OO 
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Therefore, 


A b—\ yoo 

<f>(0 = a/2 J2 h k ^( 2x - k)e' ik dx. 
k=o J ~°° 


Let y = 2x — k which implies dx = dy/2 to get 

= 4 tX> r 

k=0 J-oo 
L-l 


or simply 


1 ^ 1 , * fOO . * 

v ^ Jb=0 J -°° 


to = fnamm, 


recalling that //(£) = ^ Er=o h k e' k t. Furthermore, we get tHO = tf(£/2)tf(£/4)<£(£/4). 
This implies, 

00 £ 

m = to n « (§■ ). 


j = 1 


but <^(x) is normalized, <^(0) = / <f>(x)dx = 1. Therefore, 

00 ^ 

«« = n «(|)- 

r=i ^ 


Theorem: The Fourier Transform of {r/} is of the form, 

r( 0 = if + c£ 2M+1 + fe.o.t., 

where c G C is some constant, and ’h.o.t.’ denotes higher-order terms and will be 
used again in the proof. 

Proof: Begin with the expression for {n}: 

/ oo d 

<t>( x - y)-r-<t>(x)dx. 

-00 CLX 

If we define, 

f(x) = <f>(-x), 
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and 

9 (*) = ^(*)> 

then 

p(y) = f*9(y) 

where V is the convolution operator. The convolution theorem states that the Fourier 
Transform (continuous or discrete) of a convolution is the product of the Fourier 
Transforms: 

m = Korn- 


If we define r/ as, 


r, = p(l) 

then the semi-discrete Fourier transform of r/ is 

OO 

HO = 2Z p(£ + 2irk). 

k=—oo 


where, 

p(0 = [ p{x)e'**dx, 

J — OO 

and 

OO 

HO = H r ^ tki - 

oo 

Calculate the needed Fourier Transforms to get, 


/to = ko, 


where • denotes conjugation, and 

m = mo- 

Combine these results to get the Fourier Transform of {pi}: 

p( f) = Wf)| ! if- 
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A 

Now, we need to know the behaviour of |<^(£)|. Recall from the definitions that, 

m = (1 + e*)"(i ) u Q(e% 

where Q(e'^) does not have poles or zeros at £ = ir, see [2]. That is, //(£) has a zero 
of order M at £ = t. Therefore, H(( + 7r) has a zero of order M at £ = 0, i.e., 

H{( + w) = c( M + h.o.t., 

then 

+ 7r )| 2 = a C M + h.o.t; 

where a = |c| 2 . Recall from the definitions that, 

itf(oi 2 +itfu+*)r = i. 

Combine the two previous relations to get, 

\H(t)\ 2 = l -a? M + h.o.t. 

That is, 

^rlff«)l J le=o = 0, 

for n = 1, ...,2M — 1. The Fourier Transform of <f>{x) was found above: 

oo c 

m = n H(b 

j=i L 

We get an expression for |^(£)| 2 from, 

mW) = ft «(4) ft H(b 

i= l Z i=i z 

or 

OO f 

woi 1 = n 

j=i 1 

Now, derivatives of this expression have the form, 


and one can see that if, j^r|//(£)| 2 | 4 =o = f° r n = !>•••> 2Af — 1, then 

^won<«,=o, 

for n = 1, ..., 2M — 1. From this information we can see that a series expansion of 
l^(0| 2 & hout ( — 0 would be of the form, 

m)\* = a + K 2M + k.°.t. 

But, <f>(x) is normalized implying that <^(0) = 1 and therefore |<^(0)| 2 = 1. The 
expansion becomes, 

\m \ 2 = i + + h.o.t., 

where b € C. Recall that we are looking for the semi-discrete Fourier transform of 
r(£), which we see from above is, 

OO oo A 

r(0 = X p{(, + 27rfc) = X) iU + 2nk )\H( + 2irk )\ 2 - 

k=- oo k=-oc 

We now need to find the behaviour of + 27rfc)| 2 when k / 0. Note that in the 
expression, 

+ 2 rt )| 2 = ft I 

7=1 

the arguement, 

( + 2irk £ kir 
2 ? ~ + 2 J_1 ’ 

will for some j be equal to 

£ -I- 27 rk £ 

: — — r ~r 7T, 

2^ 2 J 

modulus 27r. That is, if & is odd in the expression (fc7r)/(2 jl ) then stop when j = 1 
since we can subtract multiples of 2ir from (kir) without changing H since H is 27 r 
periodic. If k is even, then for some j, the number (k)/( 2 J_1 ) will be odd at which 
point we again subtract some multiple of 2ir . Consequently, in the infinite product, 

W« + 2rf)| 2 =ftl«(^^)| J , 

7=1 
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there will always be a term, when k ^ 0, on the right hand side with the form, 
| + tt)| 2 . But, from above, we found that, 

\h(^ + *)\> = o(e M ) + h.o.t. 

But this implies that for k ^ 0 that the contributions to the infinite sum, 

OO 

E |^« + 2x*)| ! 

k— — QO 

are of 0(£ 2M ). That is, 

OO 

53 + 27tA:)| 2 = 1 + b£ 2M + h.o.t. 

k= — oo 

Ultimately, we need the semi-discrete Fourier transform of {r}: 

OO 

HO = * 53 (£ + + 2ttA;)| 2 , 

k=— oo 

or 

oo oo 

HO = 53 + 2ttA:)| 2 + 2wi J2 k\<j>(£ + 2irk)\ 2 . 

k=— oo k=— oo 

We already know the behaviour of the first term on the right-hand side. The second 
term on the right-hand side cannot contribute powers of £ lower than 2 M since it 
differs from the summation in the first term only by a multiple of k which does not 
allow the low power contribution when k = 0. The final step is to illustate the second 
term on the right-hand side is an odd function implying that the lowest power of ( it 
can contribute is 2 M T 1, the first odd number past 2 M. That is, define 

/({,*) = k\fa + 2*k)\\ 

and note that in the infinite summation that there is always a term with positive 
k and a term with negative k. The summation of all such +k terms and —k terms 
yields odd functions. 


m *) + no -k) = n-o k) + -o, 


38 



and this implies the desired result leaving us with the conclusion of the proof that, 

r (£) = + ib{ 2M+ ' + h.o.t. 


This completes the proof./ / 


Lemma: Let {r/} be a finite set of coefficients. These coefficients can be called 
the coefficients of a finite difference approximation to a first derivative, or these 
coefficients can be called a finite impulse response filter, or FIR filter. Furthermore, 
let the coefficients be antisymmetric: rj = — r_j which implies r 0 = 0. If the Discrete 
Fourier Transform, or DFT, of {r/} is of the form 

f(0 = i( + of 71 * 1 + h.o.t., (89) 

for some constant c € C, then the filter {r/} when applied to evenly-spaced samples 
of a function can differentiate in a finite difference sense with accuracy of order m. 
That is, {r/} can differentiate polynomials exactly up to x m . 

Before the proof, note that the DFT of a filter which is designed to approximate 
differentiation in the space domain should approximate i£ in the frequency domain: 

4~e* x = i(e iix . (90) 

ax 

That is, differentiation filters are attempting to approximate the frequency of a pure 
sinusoidal mode. 

Proof: 

Let the DFT for {r; } be defined as, 

r(0 = Yh ( 91 ) 

i 

Break up the summation to write as, 

r(0 = r 0 + Y, r < e ' V + r ‘ e ' V ■ ( 92 ) 

(>1 K-l 
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Now, impose the antisymmetry to get, 


*(()=£ n(e' fl -c *'). 


(93) 


/>1 


Using the series expansion about zero for the complex exponential one gets, 


f m-VrlfW (-«0\ 
r (f) - 2^ r K2^-p u )> 

/>1 fc =0 K - K - 


or factor to get, 


(>1 fc=0 

Interchange the summations to get, 


;s.-» l — n 


(94) 


(95) 


Jb=0 ** l> 1 


(96) 


Recall that the hypothesis was that r(^) = + c^ TO+1 + h.o.t. This implies that, 


o = !>('*- (-0‘) (97) 

for A: = 0 and for 2 < A: < m. Furthermore, 1 = X3/>i r i(l k ~ (— l) k ) must hold when 
A; = 1. But these are exactly the conditions that must hold for a filter to be able 
to differentiate exactly polynomials up through order m which are centered at zero. 
The proof for polynomials centered at some arbitrary position requires a shift in the 
index but the results are the same. This completes the proof./ / 


This completes the theorems which are at the heart of the paper. The next section 
discusses the high order of accuracy of the coefficients {r/}. 

4.4.3 ‘Superconvergence’ 

Note that the matrix RP can differentiate exactly polynomials up to degree 2 M for the 
Daubechies wavelet Dim when thought of as a finite-difference operator, even though 
the scaling function subspace Vo can only represent exactly polynomials only up to 
degree M — 1. A Similar phenomenon is encountered in the finite element method 
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under particular choices of the approximation grid and is known as superconvergence 

[ 11 ]. 

To understand the source of superconvergence in the wavelet derivative it is suf- 
ficient to have a good understanding of the proofs in the previous subsection. Let us 
note the sources of the powers of £ in the expression for the DFT of the coefficients 

{n}: 

f ( 0 = + c£ 2M+1 + h.o.t. 

Recall the definition of {r;}, 

/ GO d 

<f>( x ~ 

as well as the definitions of 4>(x) and 

H x ) = X! ^( 2 x - Oi 
/=0 

4 >{ x ) = 2 ^( 2x - 0- 

/=0 

The sources of the powers of ^ are now apparent: M powers come from <f>(x) and M + 1 
powers come from ^ <j>{x ). The ’superconvergence’ for the wavelet derivative can be 
explained by the similarity between the equations which define <f>( x) and ~t<i>{x). That 
is, they are defined by dilation equations which differ only by a multiple of 2. 

The next section will summarize and conclude this paper. 
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5 Conclusion 


A restatement of the thesis of this paper is given first followed by a brief outline of 
the argument. 

Given the evenly-spaced samples of a periodic function, /, then the matrix RP 
derived for the Daubechies wavelet D^m has the effect, when applied to /, of a finite- 
difference derivative operator of degree 2 M. 

The heart of the argument of this paper is contained in §(3) and §(4). In §(3) it 
was established that if given the evenly-spaced samples of a periodic function f(x) 
then the scaling function coefficients sq of the function at the finest scale can be 
approximated by a quadrature formula which in matrix form, 

<?o = Cf, 

yields a circulant matrix (7, where <7q approximates s"o. Furthermore, in §(3) it was 
noted that all circulant matrices with the same dimensions commute. In §(4) it was 
noted that the coefficients which map the scaling function coefficients at the finest 
scale of a periodic function to the scaling function coefficients at the finest scale of the 
derivative of the function is also circulant in form when written in matrix notation, 

So = RPs 0 . 

Furthermore, it was observed that the matrix RP can differentiate evenly-spaced sam- 
ples of a polynomial in a finite-difference sense exactly up to the order of the wavelet. 
Also, when RP is applied to the evenly-spaced samples of a periodic function then RP 
is circulant. Now, combine the results of §(3) and §(4) to get the following relation: 

/' = C~ x RPCf. 

Throughout the paper it has been noted that C and RP are circulant in form when 
f(x) is periodic and that circulant matrices commute. Therefore, the previous relation 
simply becomes, 

/' = 
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where RP is acting as a finite difference operator. 

A note concerning notation is in order. In the introduction the matrices (7, D , 
and the differentiation matrix 2? were defined. Under the scenario developed in this 
paper, the wavelet matrix C is the same as the matrix C from the introduction. The 
matrix D from the introduction becomes the matrix RP. Likewise, the matrix T> is 
also R° since for wavelets C and BP commute. That is, for evenly-spaced samples 
and periodicity B° is the wavelet differentiation matrix which has the effect of a finite 
difference operator. 

The importance of the thesis of this paper is that under periodicity and an evenly- 
spaced grid one can understand the wavelet differentiation matrix in terms of a finite 
difference operator with the accuracy given by the superconvergence theorem. 
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Appendix A: Wavelets Supported on (0,3M) 

In this appendix our wavelets are supported on [0, 3 M] where M is the number of 
vanishing moments of the wavelet. These are not the usual Daubechies wavelets, but 
for these wavelets the scaling function coefficients of a periodic function f(x) can be 
approximated with error of order M simply by sampling f(x) at the correct location. 

To begin, assume that there exist a unique tm , fixed for a fixed number of vanishing 
moments, Af, of the wavelet, such that 

J <f>(x + T\f^x m dx — 0 

for m = 1,2, — 1. Furthermore, recall the definition of the scaling function 

coefficient and expand the integrand f(x) in a Taylor series about x 0 (/q = f f (x o)): 

s° k — J f(x)<f>(x — k)dx = 

fo J <f>(x - k)dx + /o j (x - x 0 )<f>(x - k)dx + /q j (x - x 0 ) 2 <f>(x - k)dx + .... 

Now, shift the variable of integration by y = x — t — k, and choose the point of 
expansion, Xo, to be r + k to get, 



/(r + k) J <f>(y + r)dy + f'(y + k) J y<f>(y + r)dy + f"{y + k) J y 2 <f>(y + r)dy + .... 

Now, rename r as tm and the above integrals are of the form, 

J <f>(x + T M )x m dx = 0, 
and therefore vanish for m = 1, ...,M — 1 leading to, 

s° k = /(t m + k) + / (M) (t m + k) J y M <f>(y + r M )dy + 

i.e., the approximation of the scaling function coefficient s® up to order M is made 
by sampling f(x) at the position tm + k. 
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Note that all of the above calculations could have been carried out for the first 
derivative of f(x) giving an approximation to the scaling function coefficients, s'°, of 
/'(*): 

i°k = /(T + *) + /<" +,) ( t + k)J + r )iy + .... 

It was assumed above that there exist one tm such that 

J <f>(x + r M )x m dx = 0, 

for m = 1 , M — 1 . For m = 1 this Tm is easy to find: 

/ <£(* + r M )xdx = / <j>(x)(x - T M )dx 

= J x(j>(x)dx — tm J <f>(x)dx. 

But f <j>(x)dx = 1, therefore, 



That is, tm is simply the first moment of 4>(x). To find tm for m > 1 the calculations 
are simple but a bit longer and require the result from the following theorem to show 
that there is one tm which is the same for all m = 1, ..., M — 1. 

If / <f>{x)dx = 1 and there exists r such that / (f>(x+T)x m dx — 0 for m = 1, ..., M— 1 
then / <f>(x)x m dx = (/ (f>(x)xdx) m for m = 1, ..., M — 1. 

Proof: Start with 

J <j>(x + r)x m dx — 0, 

and let y = x + r to get, 

J <i>{y){y - T ) m = o. 

Using the binomial theorem this becomes, 

/ Hi i) £ ^ 7 ) y T (~ T ) m ~ Td y = °- 
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Let the moments of <j)(x) be denoted by Mi = f (f>(x)x l dx to get 


rn 


(- T ) m - r M r = 0. 


r= 0 


A simple calculation yields r = M\. Using this value of r and summing only up to 
m — 1 the previous expression becomes, 

771—1 


E 

r=0 


m 


M t + M m = 0. 


Or, 


= - E ( 7 ) (-ir~'(A ur-’Mr. 

From the hypotheses it is known that Mo = / <f>(x)dx = 1. Therefore, M p = Mf for 
p — 0, 1, and with this knowledge it is easy to show that M p = Mf for p — 2: 


m— 1 


= - E I T ] (-ir- r (M, 


r =0 


which holds for m = 1,2. Combine the powers of A/i to get, 
But, this is nothing more than, 

M m = -Mr[(i-ir-i], 


or simply, 


M m = M™, 

where m = 1,2. The proof is complete, since higher powers of m can be found by 
induction. 
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Appendix B: Moments of the Scaling Function 

In this appendix the moments of <f>(x ) will be calculated in closed form. Begin 
with the definition of the scaling function, 


= Y^^( 2x ~ k )- 

k 

Next, calculating the m-th moment of <j>(x) yields, 


J <f>(x)x m = ^hk J <f>( 2x — k)x m dx. 


Change the variable of integration such that y = 2x — k to get, 



E 


^/«y)(l/2Hi/ + *ri/2iy, 


= (i/2)“+> J <%)(»+ trd,. 

Now, recall the binomial theorem to get, 


J 4>(x)x m = (i/2) m+1 J ( 7 ) ylkm ~‘ d y 

Rewrite the moments of <f>(x) as Mi = f x , <f>(x)dx to get, 

M m = (1/2 r +1 z( 7 

1=0 V / Jt 

Now let (ii = k kk l to get 

= ( 1/2 r+' E ( 7 ) n *-< M < ■ 

Now, M m can be defined in terms of Mi for i = 0, ..., m — 1: 

M m (2" +1 - 2) = ( 7 ) *»-< M ‘ ■ 

Note that the moments /x/ can be found by direct calculation given that the Daubechies 
filter coefficients, /&*., are already known. 

The moments of <f>(x) can now be used to find the mapping from the evenly- 
spaced samples of a function f(x) to the scaling function coefficients. In section 3 
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i Mi /x, c, 

0 1 2 .3661 . 

1 .63395 1.2679 .6340 


Table 7: Scaling function moments for the Daubechies 4 wavelet. 

i Mi m Cj 

0 1 2 -.0235 

1 1.0054 2.0108 1.0426 . 

2 1.0108 2.0216 -.0199 

3 .90736 .5078 .0009 


Table 8: Scaling function moments for the Daubechies 8 wavelet. 

this mapping was denoted in matrix form as the matrix C. The elements c, which 
define the rows of this matrix have already been given for the wavelet Dg. The 
comparable results for the Z) 4 and Dg wavelets are given in the accompanying tables. 


51 







Form Approved 
OMB No 0704-0188 


Public reporting burden for this collection of information is estimated to average 1 hour per response, including the time for reviewing instructions, searching existing data sources, 
gathering and maintaining the data needed, and completing and reviewing the collection of information Send comments regarding this burden estimate or any other aspect of this 
collection of information, including suggestions for reducing this burden, to Washington Headquarters Services, Directorate for Information Operations and Reports, 1215 Jefferson 
Davis Highway, Suite 1204, Arlington, VA 22202-4302, and to the Office of Management and Budget, Paperwork Reduction Project (0704-0168), Washington, DC 20503. 


REPORT DOCUMENTATION PAGE j 


1. AGENCY USE ONLY(leave blank) I 2. REPORT DATE 

I December 1993 


4. TITLE AND SUBTITLE 

ON THE DAUBECHIES-BASED WAVELET 
DIFFERENTIATION MATRIX 


6. AUTHOR(S) 

Leland Jameson 


3. REPORT TYPE AND DATES COVERED 

Contractor Report 


5. FUNDING NUMBERS 

C NAS1-19480 
WU 505-90-52-01 


7. PERFORMING ORGANIZATION NAME(S) AND ADDRESS(ES) 

Institute for Computer Applications in Science 
and Engineering 

Mail Stop 132C, NASA Langley Research Center 
Hampton, VA 23681-0001 


8. PERFORMING ORGANIZATION 
REPORT NUMBER 

ICASE Report No. 93-95 


9. SPONSORING/MONITORING AGENCY NAME(S) AND ADDRESS(ES) 

National Aeronautics and Space Administration 
Langley Research Center 
Hampton, VA 23681-0001 


10. SPONSORING/MONITORING 
AGENCY REPORT NUMBER 

NASA CR-191583 
ICASE Report No. 93-95 


11. SUPPLEMENTARY NOTES 

Langley Technical Monitor: Michael F. Card 
Final Report 

Submitted to SIAM Journal on Scientific Computing 


12a. DISTRIBUTION/AVAILABILITY STATEMENT 

U nclassified-U nlimited 


12b. DISTRIBUTION CODE 


Subject Category 64 


13. ABSTRACT (Maximum 200 words) 

The differentiation matrix for a Daubechies- based wavelet basis will be constructed and ‘superconvergence’ will be 
proven. That is, it will be proven that under the assumption of periodic boundary conditions that the differentiation 
matrix is accurate of order 2 M } even though the approximation subspace can represent exactly only polynomials up 
to degree M — 1, where M is the number of vanishing moments of the associated wavelet. It will be illustrated that 
Daubechies- based wavelet methods are equivalent to finite difference methods with grid refinement in regions of the 
domain where small-scale structure is present. 


14. SUBJECT TERMS 

differentiation matrix; superconvergence; wavelet approximation 


IS. NUMBER OF PAGES 

53 


17. SECURITY CLASSIFICATION 
OF REPORT 

Unclassified 


NSN 7540-01-280-5500 


16. PRICE CODE 

A04 


18. SECURITY CLASSIFICATION! 19. SECURITY CLASSIFICATION 20. LIMITATION 
OF THIS PAGE I OF ABSTRACT OF ABSTRACT 

Unclassified | 


Standard Form 298(Rev. 2-89 


Prescribed by ANSI Std Z39-18 
298-102 


☆ u.S. GOVERNMENT PRINTING OFFICE: l**4 - 528-064/86105 


















National Aeronautics and 
Space Administration 

BULK RATE 

Langley Research Center 

POSTAGE & FEES PAID 

Mail Code 1 80 

NASA 

Hampton, V A 23681-00001 

Permit No. G-27 

Official Business 
Penalty for Private Use, S300 






